bulk-org
Raphael Kubler
7/19/2023
Introduction
Bulk RNA-seq analysis of cerebral organoid samples in conditions co-cultured (with/without microglia) and stimulation (ctr/LPS/IFNy).
CO vs COiMG
Organoids co-cultured with vs without microglia.
Variance analysis
PCAplot(pca_cor3_mod1, "COiMg", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='COiMG',y='Phagocytosis', fill='COiMG')) +
stat_compare_means(method = "t.test", label.y = 10, label.x = 1.3) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Phagocytosis differential PC expression",x="COiMG", y = "PC1 Phagocytosis") +
theme_classic()
Heatmap
#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.meta$COiMg[ctr.meta$COiMg %in% 'yes'] <- 'COiMg'
ctr.meta$COiMg[ctr.meta$COiMg %in% 'no'] <- 'CO'
#horizontal
ra <- rowAnnotation(
COiMG = ctr.meta$COiMg[ctr.meta$condition %in% 'ctrl'],
col = list(
COiMG = c("COiMg"='tomato','CO'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr <- t(batch.rem3_mod1[,meta3$condition %in% 'ctrl'])
hm_counts <- scale(ctr[,colnames(ctr) %in% gene_panels_subset$Microglia])
ComplexHeatmap::Heatmap(hm_counts,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(hm_counts), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
DEA
Volcano
volcano_plot <- function(res, title = NULL, subtitle = NULL, annotate_by = NULL, type ='ALS', ymax1 = 7.5, ymax2 = 8, xmax1 = -4.5, xmax2 = 4.5){
res <-
mutate(res,
sig = case_when(
padj >= 0.05 ~ "non_sig",
padj < 0.05 & abs(log2FoldChange) < 1 ~ "sig",
padj < 0.05 & abs(log2FoldChange) >= 1 ~ "sig - strong"
)) %>%
mutate(direction = ifelse(log2FoldChange > 0, "up", "down")) %>%
mutate(log2FoldChange = case_when(
log2FoldChange > 10 ~ Inf,
log2FoldChange < -10 ~ -Inf,
TRUE ~ log2FoldChange
)) %>%
mutate(class = paste(sig, direction))
if( type == "ALS"){
xpos <- 0.5
ymax <- ymax1
xlim <- c(xmax1,xmax2)
}else{
xpos <- 0.025
ymax <- ymax2
xlim <- c(-0.042, 0.042)
}
de_tally <- group_by(res, sig, direction, class) %>% tally() %>%
filter(sig != "non_sig") %>%
mutate( position = ifelse(sig == "sig", xpos, 2) ) %>%
mutate(position = ifelse( direction == "down", -1 * position, position)) %>%
mutate(n = formatC(n, format="f", big.mark=",", digits=0))
plot <- res %>%
mutate( pvalue = ifelse( pvalue < 1e-90, Inf, pvalue)) %>% #threshold at 1e16
ggplot(aes(x = log2FoldChange, y = -log10(pvalue))) +
#geom_point(aes(colour = class ), size = 0.5) +
rasterise(geom_point(aes(colour = class ), size = 0.5), dpi = 300) +
scale_colour_manual(values = c("non_sig up" = "gray",
"non_sig down" = "gray",
"sig up" = "#EB7F56",
"sig - strong up" = "#B61927",
"sig down" = "#4F8FC4",
"sig - strong down" = "dodgerblue4"
)) +
theme_bw() +
labs(y = expression(-log[10]~P~value), x = expression(log[2]~"(fold change)"), title = title, subtitle = subtitle) +
guides(colour = FALSE) +
scale_y_continuous(expand = c(0,0), limits = c(0,ymax)) +
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black")
) +
geom_text(fontface = "bold", data = de_tally, aes(x = position, y = ymax - 0.5, label = n, colour = class), size = 4) +
scale_x_continuous(limits = xlim)
if(!is.null(annotate_by)){
plot <- plot +
ggrepel::geom_text_repel(
fontface = "italic",
data = filter(res, symbol %in% annotate_by),
aes(x = log2FoldChange, y = -log10(pvalue), label = symbol),
min.segment.length = unit(0, "lines"),
size = 2.3) +
geom_point(
data = filter(res, symbol %in% annotate_by), size = 0.8, colour = "black"
) +
geom_point(aes(colour = class ),
data = filter(res, symbol %in% annotate_by), size = 0.6
)
}
return(plot)
}
res_COiMG$symbol <- rownames(res_COiMG)
volcano_plot(data.frame(res_COiMG), title = 'CO vs COiMg',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
Top 50 DEGs heatmap
top50.coimg <- rbind(results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),][results$COiMG[order(results$COiMG$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names <- rownames(top50.coimg)
ra2 <- rowAnnotation(
COiMG = meta3$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra2,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
kable(cbind('DEG'=rownames(top50.coimg),'lFC'=round(top50.coimg$log2FoldChange,2),'Cilia-associated'=rownames(top50.coimg) %in% organoid$CBC), format = "html") %>% kable_styling("striped") %>% scroll_box(width = "100%")| DEG | lFC | Cilia-associated |
|---|---|---|
| APOC2 | 8.03 | FALSE |
| GOLGA6D | 7.78 | FALSE |
| VAV1 | 7.31 | FALSE |
| FCGR3A | 7.29 | FALSE |
| C1QC | 7.14 | FALSE |
| ITGB2 | 7.11 | FALSE |
| CD86 | 7.08 | FALSE |
| OLR1 | 6.97 | FALSE |
| SPP1 | 6.86 | FALSE |
| WDFY4 | 6.6 | FALSE |
| ALOX5AP | 6.58 | FALSE |
| TYROBP | 6.55 | FALSE |
| ADAMDEC1 | 6.55 | FALSE |
| BIN2 | 6.54 | FALSE |
| CD53 | 6.49 | FALSE |
| MS4A7 | 6.45 | FALSE |
| SRGN | 6.4 | FALSE |
| CCRL2 | 6.36 | FALSE |
| CD52 | 6.35 | FALSE |
| FCER1G | 6.29 | FALSE |
| LSP1 | 6.27 | FALSE |
| PECAM1 | 6.25 | FALSE |
| CD48 | 6.19 | FALSE |
| DCSTAMP | 6.19 | FALSE |
| CHI3L1 | 6.14 | FALSE |
| CD164L2 | -4.48 | FALSE |
| STOML3 | -4.39 | TRUE |
| MMP10 | -4.36 | FALSE |
| GIPC2 | -4.21 | FALSE |
| LRRC18 | -4.17 | FALSE |
| FAM166C | -4.14 | FALSE |
| LRRC71 | -3.96 | FALSE |
| IL5RA | -3.7 | FALSE |
| TTC22 | -3.67 | FALSE |
| CD302 | -3.66 | FALSE |
| C20orf85 | -3.66 | FALSE |
| CLDN2 | -3.6 | FALSE |
| MFRP | -3.56 | FALSE |
| TACR3 | -3.42 | FALSE |
| CIBAR2 | -3.42 | FALSE |
| C11orf97 | -3.41 | FALSE |
| F5 | -3.38 | FALSE |
| SAG | -3.34 | FALSE |
| ADGB | -3.31 | FALSE |
| TEKT4 | -3.29 | FALSE |
| VWA5B1 | -3.28 | FALSE |
| DEUP1 | -3.28 | FALSE |
| DYNLT5 | -3.21 | FALSE |
| ARSI | -3.2 | FALSE |
| FAM81B | -3.18 | TRUE |
GO
dotplot(all_res$COiMG$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$COiMG$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$COiMG$DEGs[,1]),"Down-regulated"=na.omit(all_res$COiMG$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
cgenes <- readxl::read_xlsx('Gene lists for organoid paper.xlsx', col_names = T)
m <- data.frame(t(toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`))[toupper(na.omit(cgenes$`ADHD GWAS associated genes through FUMA`)) %in% rownames(bulkorg_counts3)]), check.names = F)
for(i in names(cgenes)[-1]){
k <- na.omit(cgenes[i])
s <- k[as.matrix(k)[,1] %in% rownames(bulkorg_counts3),]
m <- rbind.fill(data.frame(m),data.frame(t(s)))
}
rownames(m) <- names(cgenes)
m <- t(m)
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
IFNy stimulation
Let’s do this for IFNy stimulation
Variance analysis
PCAplot(pca_cor3_mod1, "condition", Shape = "Line", PoV.df=PoV_cor3_mod1, pc.1 = 1, pc.2 = 2) #M30 & M33 are outleirs, let's remove them for now
ggplot(pca.df_upgraded,aes_string(x='Condition',y='`IFNy`', fill='Condition')) +
#stat_compare_means(method = "t.test", label.y = 20, label.x = 2) +
geom_boxplot() +
geom_boxplot(width=0.1, fill="white") +
scale_fill_brewer(palette="Dark2") +
labs(title="Condition differential PC expression",x="Condition", y = "PC1 IFNy response") +
theme_classic()
DEA
Volcano
res_IFNg$symbol <- rownames(res_IFNg)
volcano_plot(data.frame(res_IFNg), title = 'CTR vs IFNy',
annotate_by=c('CD36','CD68','CX3CR1','CSF1R','TREM2','TMEM119','ITGAX'), ymax1 = 30, ymax2 = 31)
COiMg.degs <- all_res$COiMG$DEGsTop 50 DEGs heatmap
top50.ifng <- rbind(results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = T),]$padj < 0.05,][1:25,],
results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),][results$IFNg[order(results$IFNg$log2FoldChange, decreasing = F),]$padj < 0.05,][1:25,])
top50.names2 <- rownames(top50.ifng)
ra3 <- rowAnnotation(
Stimulation = meta3$condition,
col = list(
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(scale(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra3,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% top50.names2,])), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "Z-score expression")
Which genes of the top DEGs are cilia-associated (CBC)?
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
kable(cbind('DEG'=rownames(top50.ifng),'lFC'=round(top50.ifng$log2FoldChange,2),'Cilia-associated'=rownames(top50.ifng) %in% organoid$CBC), format = "html") %>% kable_styling("striped") %>% scroll_box(width = "100%")| DEG | lFC | Cilia-associated |
|---|---|---|
| IDO1 | 9.96 | FALSE |
| GBP5 | 9.22 | FALSE |
| ETV7 | 9.05 | FALSE |
| CXCL10 | 8.94 | FALSE |
| SIGLEC7 | 8.91 | FALSE |
| GBP1 | 8.72 | FALSE |
| CXCL9 | 8.71 | FALSE |
| APOL6 | 8.63 | FALSE |
| APOL1 | 8.54 | FALSE |
| PLAAT4 | 8.34 | FALSE |
| GBP3 | 8.32 | FALSE |
| SAMD9L | 8.26 | FALSE |
| GBP2 | 8.22 | FALSE |
| XAF1 | 8.07 | FALSE |
| PSMB9 | 8.03 | FALSE |
| CXCL11 | 8.02 | FALSE |
| GBP6 | 7.77 | FALSE |
| CD274 | 7.57 | FALSE |
| OAS2 | 7.56 | FALSE |
| DES | 7.5 | FALSE |
| UBD | 7.49 | FALSE |
| OR2I1P | 7.45 | FALSE |
| OAS1 | 7.43 | FALSE |
| NLRC5 | 7.41 | FALSE |
| BATF2 | 7.31 | FALSE |
| CROCC2 | -5.63 | FALSE |
| F5 | -5.55 | FALSE |
| STMND1 | -4.97 | FALSE |
| MUSK | -4.91 | FALSE |
| ARHGAP36 | -4.65 | FALSE |
| PTPRC | -4.59 | TRUE |
| IGF1 | -4.58 | FALSE |
| KL | -4.54 | FALSE |
| STOML3 | -4.52 | TRUE |
| ATP13A5 | -4.47 | FALSE |
| SLC39A12 | -4.4 | TRUE |
| CYTL1 | -4.35 | FALSE |
| C11orf97 | -4.26 | FALSE |
| SERPINA3 | -4.15 | FALSE |
| TMEM125 | -4.13 | FALSE |
| FAM81B | -4.08 | TRUE |
| C20orf85 | -4.03 | FALSE |
| CAPSL | -3.99 | TRUE |
| CIBAR2 | -3.96 | FALSE |
| TRPM1 | -3.95 | FALSE |
| C10orf105 | -3.92 | FALSE |
| HMGCS2 | -3.87 | FALSE |
| CFAP73 | -3.81 | FALSE |
| PMCH | -3.75 | FALSE |
| MMP10 | -3.69 | FALSE |
GO
dotplot(all_res$IFNg$GO.up, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
dotplot(all_res$IFNg$GO.down, split="ONTOLOGY", showCategory=10) + facet_grid(ONTOLOGY~., scale="free")
Enrichment
For organoid cell-types.
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
ogenes <- readxl::read_xlsx('gene list_cell type.xlsx', col_names = T)
organoid <- ogenes[,colnames(ogenes) %in% c('CN 1','CN 2','CN 3','CN 4','CN 5','IN','Neuron','Inter','OPC/OL','GPC','NEC 1',
'NEC 2','NEC 3','NEC 4','AS1','AS2','AS3','PGC 1','PGC 2','BRC','CBC','UPRC1','UPRC2')][-1,]
#gene set enrichment analysis:
DEG.enrich.res <- matrix(nrow=2*length(colnames(organoid)), ncol=6)
DEG.rep <- NULL
for (i in names(organoid))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, organoid[,i], universe=35324))
for (x in colnames(organoid)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(organoid)[,x]), universe=na.omit(rownames(res))))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(organoid)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
For NDD gene lists.
setwd('/Users/kubler01/Documents/R projects/gene lists/ndd')
#gene set enrichment analysis:
ndd.enrich <- matrix(nrow=2*length(colnames(m)), ncol=6)
ndd.rep <- NULL
for (i in colnames(m))
ndd.rep <- c(ndd.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
rownames(ndd.enrich) <- ndd.rep
colnames(ndd.enrich) <- colnames(GSEA.byMod(mod.gl=DEGs, m[,1], universe=35324))
for (x in colnames(m)){
ndd.enrich[gsub("*_up","",gsub("*_down","",rownames(ndd.enrich))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(as.matrix(m)[,x]), universe=na.omit(rownames(res))))
}
ndd.enrich <- data.frame(ndd.enrich)
ndd.enrich$log.q <- -log2(ndd.enrich[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
ndd.enrich.ress <- cbind("Upregulated"=ndd.enrich[c(seq(1, length(rownames(ndd.enrich))-1, by = 2)),][,7],
"Downregulated"=ndd.enrich[c(seq(2, length(rownames(ndd.enrich)), by = 2)),][,7])
rownames(ndd.enrich.ress) <- colnames(m)
ComplexHeatmap::Heatmap(ndd.enrich.ress,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
DEGs$`Down-regulated`[DEGs$`Down-regulated` %in% as.character(na.omit(m[,12]))]## [1] "CYP7B1" "GRIN2A"For microglia pathways.
DEG.enrich.res <- matrix(nrow=2*length(colnames(gene_panels_subset)), ncol=6)
DEG.rep <- NULL
for (i in names(gene_panels_subset))
DEG.rep <- c(DEG.rep, paste(i, "_up", sep=""), paste(i, "_down", sep=""))
DEGs <- list("Up-regulated"=na.omit(all_res$IFNg$DEGs[,1]),"Down-regulated"=na.omit(all_res$IFNg$DEGs[,2]))
rownames(DEG.enrich.res) <- DEG.rep
colnames(DEG.enrich.res ) <- colnames(GSEA.byMod(mod.gl=DEGs, gene_panels_subset[,i], universe=35324))
for (x in colnames(gene_panels_subset)){
DEG.enrich.res[gsub("*_up","",gsub("*_down","",rownames(DEG.enrich.res))) %in% x,] <- as.matrix(GSEA.byMod(mod.gl=DEGs, na.omit(gene_panels_subset[,x]), universe=rownames(res)))
}
DEG.enrich.res <- data.frame(DEG.enrich.res)
DEG.enrich.res$log.q <- -log2(DEG.enrich.res[,6])
#DEG.enrich.res[is.infinite(DEG.enrich.res$log.q),]$log.q <- 1
DEG.enrich <- cbind("Upregulated"=DEG.enrich.res[c(seq(1, length(rownames(DEG.enrich.res))-1, by = 2)),][,7],
"Downregulated"=DEG.enrich.res[c(seq(2, length(rownames(DEG.enrich.res)), by = 2)),][,7])
rownames(DEG.enrich) <- names(gene_panels_subset)
DEG.enrich[DEG.enrich > 20] <- 20
ComplexHeatmap::Heatmap(DEG.enrich,
cluster_rows = F,
cluster_columns = F,
show_row_dend = F,
show_column_dend = F,
show_row_names = T,
col = colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"))(30),
name = "log2(q)",
)
Cell-type marker analysis
Cell-types with only 1 gene found back in our expression data carry principal component expression = 0.
Heatmaps of gene & PC expression
setwd('/Users/kubler01/Documents/R projects/gene lists/organoids')
cell.types <- as.matrix(readxl::read_xlsx('Typical markers cell types organoids 07182023.xlsx', col_names = F))
PC.covs <- matrix(ncol=length(cell.types[,1]),nrow=length(colnames(batch.rem3_mod1)))
colnames(PC.covs) <- cell.types[,1]
for(i in cell.types[,1]){
pca.cov <- prcomp(t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],]))
PC.covs[,i] <- pca.cov$x[,1]
#if (table(rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1])[[2]] == 1){
# PC.covs[,i] <- t(batch.rem3_mod1[rownames(batch.rem3_mod1) %in% cell.types[cell.types[,1] %in% i,][-1],])
#}
}
rownames(PC.covs) <- colnames(batch.rem3_mod1)
ra4 <- rowAnnotation(
Stimulation = meta3$condition,
COiMG = meta3$COiMg,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ComplexHeatmap::Heatmap(PC.covs,
cluster_rows = F,
cluster_columns = F,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ra4,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(PC.covs), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's do that for ctr/ifng and coimg/CO only
#for control only
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
#horizontal
ctr.ra <- rowAnnotation(
COiMG = ctr.meta$COiMg,
col = list(
COiMG = c("yes"='tomato','no'='lightblue')),
show_annotation_name = T,
show_legend = T)
ctr.pca <- PC.covs[meta3$condition %in% 'ctrl',]
ComplexHeatmap::Heatmap(ctr.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
ifng.meta <- meta3[!meta3$condition %in% 'LPS',]
#horizontal
ifng.ra <- rowAnnotation(
COiMG = ifng.meta$COiMg,
Stimulation = ifng.meta$condition,
col = list(COiMG = c('yes'='black','no'='white'),
Stimulation = c("ctrl"='lightblue','IFNg'='tomato', 'LPS'='darkred')),
show_annotation_name = T,
show_legend = T)
ifng.pca <- PC.covs[!meta3$condition %in% 'LPS',]
ComplexHeatmap::Heatmap(ifng.pca,
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ifng.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(ctr.pca), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
#Let's check only choroid plexus & ependymal
cp.ep <- c(na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][1,][-1]),
na.omit(cell.types[cell.types[,1] %in% c('Chorioid plexus','Ependymal'),][2,][-1]))
ctr.df <- batch.rem3_mod1[,colnames(batch.rem3_mod1) %in% rownames(ctr.meta)]
ComplexHeatmap::Heatmap(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,])),
cluster_rows = T,
cluster_columns = T,
show_row_dend = T,
show_column_dend = T,
show_row_names = F,
show_column_names = F,
left_annotation = ctr.ra,
bottom_annotation = HeatmapAnnotation(
text = anno_text(colnames(scale(t(ctr.df[rownames(ctr.df) %in% cp.ep,]))), rot = 300, just = "left",gp=gpar(fontsize=8))),
col = colorRampPalette(rev(RColorBrewer::brewer.pal(9,"RdBu")))(10),
name = "PC1 expression")
Boxplots
Let’s make boxplots instead
ifng.pca <- PC.covs[!meta3$condition %in% 'LPS',]
ifng.pca_long <- reshape::melt(ifng.pca)
colnames(ifng.pca_long) <- c('Sample','Cell.type','PC')
ifng.pca_long$Stimulation <- ifng.meta$condition
ggplot(ifng.pca_long,aes_string(x='Cell.type',y='PC', fill='Stimulation')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
ctr.meta <- meta3[meta3$condition %in% 'ctrl',]
ctr.pca <- PC.covs[meta3$condition %in% 'ctrl',]
coimg.pca_long <- reshape::melt(ctr.pca)
colnames(coimg.pca_long) <- c('Sample','Cell.type','PC')
coimg.pca_long$Coculture <- ctr.meta$COiMg
ggplot(coimg.pca_long,aes_string(x='Cell.type',y='PC', fill='Coculture')) +
#stat_compare_means(method = "t.test", label.y = 6) +
geom_boxplot() +
labs(title="Cell-type marker expression",x="Cell-type", y = "PC1 marker expression") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5))